# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to run multilevel models and create:
#   - Table A.5 (multilevel of main results)
#   - Table A.13 (multilevel of Freedom House interaction)
#   - Table A.17 (multilevel of KOF interaction)

# NOTE: These models take a long time to run. 

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data
myData <- readRDS("data/afrobarometer.rds")

# Transform variables:
# all to numeric b/c simcf doesn't like factors
myData$sexuality2 <- as.numeric(myData$sexuality2)
myData$tolerance_noLGBT2 <- as.numeric(myData$tolerance_noLGBT2)
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)
myData$tolerance_noLGBT2 <- as.numeric(myData$tolerance_noLGBT2)
myData$female <- as.numeric(myData$female)
myData$education <- as.numeric(myData$education)
myData$religiosity <- as.numeric(myData$religiosity)
myData$age <- as.character(myData$age) #age needs to go through character first
myData$age <- as.numeric(myData$age) #then to numeric to maintain real values
myData$income_water <- as.numeric(myData$income_water)
myData$urban <- as.numeric(myData$urban)


# MODELS

##################
# Media Aggregate
##################


# B. media aggregate: VSVI Country + VI District
model.b.aggregate <-
  sexuality2 ~
  media_aggregate +
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + media_aggregate | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.aggregate, myData, na.rm = TRUE) #only extract data used in model
aggregate.b.result <- glmer(model.b.aggregate, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(aggregate.b.result, file = project_path("analysis/results", "aggregate.b.result.rds"))

# B. media aggregate: VSVI Country + VI District w/o tolerance
model.b.aggregate <-
   sexuality2 ~
   media_aggregate +
   female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
   (1 + media_aggregate | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.aggregate, myData, na.rm = TRUE) #only extract data used in model
aggregate.b.result.notol <- glmer(model.b.aggregate, family=binomial, data=mdata,
                             nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(aggregate.b.result.notol, file = project_path("analysis/results", "aggregate.b.result.notol.rds"))

##########
# Radio
#########

# B. radio: VSVI Country + VI District
model.b.radio <-
  sexuality2 ~
  radio +
  media_noradio + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + radio | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.radio, myData, na.rm = TRUE) #only extract data used in model
radio.b.result <- glmer(model.b.radio, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(radio.b.result, file = project_path("analysis/results", "radio.b.result.rds"))

# B. radio: VSVI Country + VI District w/o tolerance
model.b.radio <-
  sexuality2 ~
  radio +
  media_noradio + #consumption of other mediums
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + radio | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.radio, myData, na.rm = TRUE) #only extract data used in model
radio.b.result.notol <- glmer(model.b.radio, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(radio.b.result.notol, file = project_path("analysis/results", "radio.b.result.notol.rds"))

##########
# TV
##########

# B. tv: VSVI Country + VI District
model.b.tv <-
  sexuality2 ~
  tv +
  media_notv + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + tv | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.tv, myData, na.rm = TRUE) #only extract data used in model
tv.b.result <- glmer(model.b.tv, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(tv.b.result, file = project_path("analysis/results", "tv.b.result.rds"))

# B. tv: VSVI Country + VI District w/o tolerance
model.b.tv <-
  sexuality2 ~
  tv +
  media_notv + #consumption of other mediums
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + tv | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.tv, myData, na.rm = TRUE) #only extract data used in model
tv.b.result.notol <- glmer(model.b.tv, family=binomial, data=mdata,
                     nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(tv.b.result.notol, file = project_path("analysis/results", "tv.b.result.notol.rds"))

############
# Newspaper
############

# B. newspaper: VSVI Country + VI District
model.b.newspaper <-
  sexuality2 ~
  newspaper +
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + newspaper | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.newspaper, myData, na.rm = TRUE) #only extract data used in model
newspaper.b.result <- glmer(model.b.newspaper, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(newspaper.b.result, file = project_path("analysis/results", "newspaper.b.result.rds"))

# B. newspaper: VSVI Country + VI District w/o tolerance
model.b.newspaper <-
  sexuality2 ~
  newspaper +
  media_nopaper + #consumption of other mediums
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + newspaper | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.newspaper, myData, na.rm = TRUE) #only extract data used in model
newspaper.b.result.notol <- glmer(model.b.newspaper, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(newspaper.b.result.notol, file = project_path("analysis/results", "newspaper.b.result.notol.rds"))

###########
# Internet
###########

# B. internet: VSVI Country + VI District
model.b.internet <-
  sexuality2 ~
  internet +
  media_nointernet + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + internet | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.internet, myData, na.rm = TRUE) #only extract data used in model
internet.b.result <- glmer(model.b.internet, family=binomial, data=mdata, 
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(internet.b.result, file = project_path("analysis/results", "internet.b.result.rds"))

# B. internet: VSVI Country + VI District w/o tolerance
model.b.internet <-
  sexuality2 ~
  internet +
  media_nointernet + #consumption of other mediums
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + internet | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.internet, myData, na.rm = TRUE) #only extract data used in model
internet.b.result.notol <- glmer(model.b.internet, family=binomial, data=mdata, 
                           nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(internet.b.result.notol, file = project_path("analysis/results", "internet.b.result.notol.rds"))

###############
# Social Media
###############

# B. social media: VSVI Country + VI District
model.b.social_media <-
  sexuality2 ~
  social_media +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + social_media | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.social_media, myData, na.rm = TRUE) #only extract data used in model
social_media.b.result <- glmer(model.b.social_media, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(social_media.b.result, file = project_path("analysis/results", "social_media.b.result.rds"))

# B. social media: VSVI Country + VI District w/o tolerance
model.b.social_media <-
  sexuality2 ~
  social_media +
  media_nosocialmedia + #consumption of other mediums
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + social_media | country) + (1 | district)  #add VSVI for country & VI for District
mdata <- extractdata(model.b.social_media, myData, na.rm = TRUE) #only extract data used in model
social_media.b.result.notol <- glmer(model.b.social_media, family=binomial, data=mdata,
                               nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(social_media.b.result.notol, file = project_path("analysis/results", "social_media.b.result.notol.rds"))

# load all results to print
aggregate.b.result <- readRDS("analysis/results/aggregate.b.result.rds")
radio.b.result <- readRDS("analysis/results/radio.b.result.rds")
tv.b.result <- readRDS("analysis/results/tv.b.result.rds")
newspaper.b.result <- readRDS("analysis/results/newspaper.b.result.rds")
internet.b.result <- readRDS("analysis/results/internet.b.result.rds")
social_media.b.result <- readRDS("analysis/results/social_media.b.result.rds")

# Generate Latex table: A.5 in Appendix 
stargazer(aggregate.b.result, radio.b.result, tv.b.result, newspaper.b.result, internet.b.result, social_media.b.result,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on LGBT Attitudes (Multilevel Models)",
          omit= "country",
          label = "main_multilevel_table",
          notes = c("")
)

########################
# Free Press Interaction
# w/ VSVI by country
# & VI by district 
########################

#convert FH Score to numeric
#this is original score 0-100 where 100 = not free
myData$fh_score <- as.numeric(myData$fh_score) 
myData$fh_net_score <- as.numeric(myData$fh_net_score)

# D. aggregate
model.d.aggregate <-
  sexuality2 ~
  media_aggregate + media_aggregate:fh_scale +
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + media_aggregate | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.aggregate, myData, na.rm = TRUE) #only extract data used in model
aggregate.d.result <- glmer(model.d.aggregate, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit

# D. radio
model.d.radio <-
  sexuality2 ~
  radio + radio:fh_scale  +
  media_noradio+ #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + radio | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.radio, myData, na.rm = TRUE) #only extract data used in model
radio.d.result <- glmer(model.d.radio, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit

# D. tv
model.d.tv <-
  sexuality2 ~
  tv + tv:fh_scale + 
  media_notv + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + tv | country) + (1 | district) #add VSVI for country & VI for district 
mdata <- extractdata(model.d.tv, myData, na.rm = TRUE) #only extract data used in model
tv.d.result <- glmer(model.d.tv, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit

# newspaper
model.d.newspaper <-
  sexuality2 ~
  newspaper + newspaper:fh_scale + 
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + newspaper | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.newspaper, myData, na.rm = TRUE) #only extract data used in model
newspaper.d.result <- glmer(model.d.newspaper, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit

# internet
model.d.internet <-
  sexuality2 ~
  internet + internet:fh_scale +
  media_nointernet + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + internet | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.internet, myData, na.rm = TRUE) #only extract data used in model
internet.d.result <- glmer(model.d.internet, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit

# social_media
model.d.social_media <-
  sexuality2 ~
  social_media + social_media:fh_scale +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + social_media | country) + (1 | district) #add VSVI for country & VI for district 
mdata <- extractdata(model.d.social_media, myData, na.rm = TRUE) #only extract data used in model
social_media.d.result <- glmer(model.d.social_media, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit

# Generate Latex Table: A.13 in Appendix 
stargazer(aggregate.d.result, radio.d.result, tv.d.result, newspaper.d.result, internet.d.result, social_media.d.result,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          #column.labels=c("Logit", ""), # "Pol Herf"),
          title = "Effect of Media Consumption on LGBT Attitudes (Multilevel Models)",
          omit= "country",
          notes = c("")
)

########################
# Globalization Interaction
# w/ VSVI by country
# & VI by district 
########################

# D. aggregate
model.d.aggregate <-
  sexuality2 ~
  media_aggregate + media_aggregate:KOFSoGI +
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + media_aggregate | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.aggregate, myData, na.rm = TRUE) #only extract data used in model
aggregate.d.result <- glmer(model.d.aggregate, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit

saveRDS(aggregate.d.result, file = project_path("analysis/results", "aggregate.kof_interaction.result.rds"))

# D. radio
model.d.radio <-
  sexuality2 ~
  radio + radio:KOFSoGI  +
  media_noradio+ #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + radio | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.radio, myData, na.rm = TRUE) #only extract data used in model
radio.d.result <- glmer(model.d.radio, family=binomial, data=mdata,
                        nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(radio.d.result, file = project_path("analysis/results", "radio.kof_interaction.result.rds"))

# D. tv
model.d.tv <-
  sexuality2 ~
  tv + tv:KOFSoGI + 
  media_notv + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + tv | country) + (1 | district) #add VSVI for country & VI for district 
mdata <- extractdata(model.d.tv, myData, na.rm = TRUE) #only extract data used in model
tv.d.result <- glmer(model.d.tv, family=binomial, data=mdata,
                     nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(tv.d.result, file = project_path("analysis/results", "tv.kof_interaction.result.rds"))

# newspaper
model.d.newspaper <-
  sexuality2 ~
  newspaper + newspaper:KOFSoGI + 
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + newspaper | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.newspaper, myData, na.rm = TRUE) #only extract data used in model
newspaper.d.result <- glmer(model.d.newspaper, family=binomial, data=mdata,
                            nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(newspaper.d.result, file = project_path("analysis/results", "newspaper.kof_interaction.result.rds"))

# internet
model.d.internet <-
  sexuality2 ~
  internet + internet:KOFSoGI +
  media_nointernet + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + internet | country) + (1 | district) #add VSVI for country & VI for district
mdata <- extractdata(model.d.internet, myData, na.rm = TRUE) #only extract data used in model
internet.d.result <- glmer(model.d.internet, family=binomial, data=mdata,
                           nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(internet.d.result, file = project_path("analysis/results", "internet.kof_interaction.result.rds"))

# social_media
model.d.social_media <-
  sexuality2 ~
  social_media + social_media:KOFSoGI +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 + social_media | country) + (1 | district) #add VSVI for country & VI for district 
mdata <- extractdata(model.d.social_media, myData, na.rm = TRUE) #only extract data used in model
social_media.d.result <- glmer(model.d.social_media, family=binomial, data=mdata,
                               nAGQ = 0 ) # optimizer for faster but less precise fit
saveRDS(social_media.d.result, file = project_path("analysis/results", "social_media.kof_interaction.result.rds"))

# Generate Latex Table: A.17 in Appendix
stargazer(aggregate.d.result, radio.d.result, tv.d.result, newspaper.d.result, internet.d.result, social_media.d.result,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          #column.labels=c("Logit", ""), # "Pol Herf"),
          title = "Effect of Media Consumption on LGBT Attitudes (Multilevel Models)",
          omit= "country",
          label = "multilevel_KOF_table",
          notes = c("")
)